{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "29a77281-d15f-4459-b42b-90639d9fd158",
   "metadata": {},
   "source": [
    "# Using the ASE calculators\n",
    "\n",
    "This tutorial demonstrates how to use the [ASE](https://wiki.fysik.dtu.dk/ase/) calculators for neuroevolution potentials (NEP), which allow one to calculate energies, forces and stresses for an atomic configuration, specified in the form of an [Atoms object](https://wiki.fysik.dtu.dk/ase/ase/atoms.html).\n",
    "This enables one to programmatically calculate properties that would otherwise require writing a large number of GPUMD input files and potentially tedious extraction of the results.\n",
    "\n",
    "`calorine` provides two different ASE calculators for NEP calculations, one that uses the GPU implementation and one that uses the CPU implementation of NEP.\n",
    "For smaller calculations the CPU calculators is usually more performant.\n",
    "For very large simulations and for comparison the GPU calculator can be useful as well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8f7acff",
   "metadata": {},
   "source": [
    "All models and structures required for running this and the other tutorial notebooks can be obtained from [Zenodo](https://zenodo.org/record/10658778).\n",
    "The files are also available in the `tutorials/` folder in the [GitLab repository](https://gitlab.com/materials-modeling/calorine/-/tree/master/tutorials)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e69f41d6",
   "metadata": {
    "tags": []
   },
   "source": [
    "## CPU-based calculator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fb133d3",
   "metadata": {},
   "source": [
    "### Basic usage\n",
    "\n",
    "First we define an atomic structure in the form of an [Atoms object](https://wiki.fysik.dtu.dk/ase/ase/atoms.html).\n",
    "Next we create a calculator instance by specifying the path to a NEP model file in [nep.txt format](https://gpumd.org/nep/output_files/nep_txt.html).\n",
    "In this tutorial, we use a [NEP4 model for PbTe](https://github.com/brucefan1983/GPUMD/blob/ad07f7a031c87a9f235e54035d179a118ab2ff9a/examples/11_NEP_potential_PbTe/nep.txt). This model is included in the `calorine/tutorials/` folder.\n",
    "Finally, we attach the calculator to the `Atoms` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ccb1d625",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.io import read\n",
    "from ase.build import bulk\n",
    "from calorine.calculators import CPUNEP\n",
    "\n",
    "structure = bulk('PbTe', crystalstructure='rocksalt', a=6.7)\n",
    "calc = CPUNEP('nep-PbTe.txt')\n",
    "structure.calc = calc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "506e961e",
   "metadata": {},
   "source": [
    "Now, we can readily calculate energies, forces and stresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "12dbc5b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy (eV): -7.663163505693399\n",
      "Forces (eV/Å):\n",
      " [[ 9.43689571e-16 -2.08961602e-16 -4.73106562e-16]\n",
      " [-9.50628465e-16  2.08961602e-16  4.73106562e-16]]\n",
      "Stress (GPa):\n",
      " [ 1.09085876e-02  1.09085876e-02  1.09085876e-02 -2.53879747e-18\n",
      "  1.86911426e-18  8.83104401e-18]\n"
     ]
    }
   ],
   "source": [
    "print('Energy (eV):', structure.get_potential_energy())\n",
    "print('Forces (eV/Å):\\n', structure.get_forces())\n",
    "print('Stress (GPa):\\n', structure.get_stress())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81f94877",
   "metadata": {},
   "source": [
    "### Calculate energy-volume curve\n",
    "\n",
    "To demonstrate the capabilities of the ASE calculator, we will now calculate an energy-volume curve with the PbTe potential."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "861b6acb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGxCAYAAACOSdkqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABc0UlEQVR4nO3dd1zV9f4H8Nc5jMM+AoIcFJDhQhygkpJaiiZq4MgyV66GdsvymmmpmZlp3va4thwZZlmaq0JTy9w4QEUEFRFkK0eGbM75/v7gx7kSw3PgHL4czuv5eJxH8h3w5htyXn6mRBAEAUREREQmSCp2AURERERiYRAiIiIik8UgRERERCaLQYiIiIhMFoMQERERmSwGISIiIjJZDEJERERkshiEiIiIyGSZi11AS6dWq5GRkQF7e3tIJBKxyyEiIiItCIKAwsJCuLu7Qyqtv92HQeg+MjIy4OHhIXYZRERE1Ag3b95Ehw4d6j3PIHQf9vb2AKoepIODg8jVEBERkTYKCgrg4eGheR+vD4PQfVR3hzk4ODAIERERGZn7DWvhYGkiIiIyWQxCREREZLIYhIiIiMhkMQgRERGRyWIQIiIiIpPFIEREREQmi0GIiIiITBaDEBEREZksLqhIREREzU6lFhCdrEROYSlc7a0Q7O0EM2nz7+nJIERERETNKiouEyv2xCMzv1RzTCG3wvJwf4QFKJq1FqPpGouIiICnpyesrKygUCgwbdo0ZGRkNHjPc889B19fX1hbW8PFxQVjxoxBQkJCM1VMRERE/xQVl4m5kedqhCAAyMovxdzIc4iKy2zWeowmCA0ZMgTbtm1DYmIitm/fjqSkJEyYMKHBe/r06YONGzfi8uXL2LdvHwRBwCOPPAKVStVMVRMREVE1lVrAij3xEOo4V31sxZ54qNR1XWEYEkEQmu+r6dHu3bsxduxYlJWVwcLCQqt7Lly4gF69euHatWvw9fXV6p6CggLI5XLk5+dz01UiIqImOJGUi0lfn7zvdVuf6Y8Bvs5N+lravn8b5RghpVKJLVu2ICQkROsQVFRUhI0bN8Lb2xseHh71XldWVoaysjLNxwUFBU2ul4iIiICcwtL7X6TDdfpgNF1jALBo0SLY2trC2dkZqamp2LVr133v+e9//ws7OzvY2dnh999/xx9//AFLS8t6r1+9ejXkcrnm1VBoIiIiIu252lvp9Tp9EDUILV68GBKJpMHXvYObFy5ciJiYGOzfvx9mZmZ46qmncL+evSlTpiAmJgaHDx9G586d8cQTT6C0tP6k+dprryE/P1/zunnzpt6+XyIiIlMW7O0EZ9v6GyMkqJo9Fuzt1Gw1iTpG6NatW8jNzW3wGh8fnzpbcNLS0uDh4YHjx49jwIABWn298vJyODo64ptvvsGkSZO0uodjhIiIiPTjQloenvjyBEor1LXOVa8gtG5qkF6m0BvFGCEXFxe4uLg06l61uuoh3jue534EQYAgCDrdQ0RERE13MS0fU785hdIKNXxdbHG3rBLZBf97P3YTaR0hoxgsferUKZw+fRoDBw6Eo6MjkpKSsGzZMvj6+mpag9LT0xEaGorNmzcjODgY169fx48//ohHHnkELi4uSEtLw5o1a2BtbY1Ro0aJ/B0RERGZjotp+ZjyzUkUlFair5cjNs0KhrWFGVeW1paNjQ127NiB5cuXo6ioCAqFAmFhYVi6dClkMhkAoKKiAomJiSguLgYAWFlZ4ciRI/joo49w584dtGvXDoMHD8bx48fh6uoq5rdDRERkMuLS8zF1/akaIchOVhU/mjpFXh+Mdh2h5sIxQkRERI0Tl56PKd+cQn5JBfp4OeLbe0KQoWn7/m1U0+eJiIjIONwbgoI822DTzH7NFoJ0wSBEREREevXPEPTtrGDYW2m3AHJzYxAiIiIivbmUUTUmKL+kAoEtPAQBDEJERESkJ5cyqlqC8oor0Nuj5YcggEGIiIiI9CA+o6BGCNo8OxgOLTwEAQxCRERE1ERVIegk8oor0MuIQhDAIERERERNUB2C7hRXoFcHOTbPMp4QBDAIERERUSNdzvxHCJr9AOTWxhOCAAYhIiIiaoSErKoxQcYcggAGISIiItJRQlYBJn99CsqicvQ04hAEMAgRERGRDhKzCjUhqEd7Ob6bZbwhCGAQIiIiIi1VhaCTmhAUOfsByG2MNwQBDEJERESkhSvZVSEot6gcAe0dWkUIAoCWt/sZERERiUalFhCdrEROYSlc7a0Q7O2EpFt3MemrqhDU3b31hCCAQYiIiIj+X1RcJlbsiUdmfqnmWFs7S5RVqlFYWonu7g7Y8vQDaGNjKWKV+sUgRERERIiKy8TcyHMQ/nH89t1yAECHNtatLgQBHCNERERk8lRqASv2xNcKQfeqUKtb/AaqjcEgREREZOKik5U1usPqkl1QhuhkZTNV1HwYhIiIiExcTmHDIUjX64wJgxAREZGJc7W30ut1xoRBiIiIyMQFezvBzaH+kCMBoJBXTaVvbRiEiIiITJyZVIJeHvI6z0n+/7/Lw/1hJpXUeY0xYxAiIiIycdtO38S+S9kAUGvfMDe5FdZNDUJYgEKM0gyO6wgRERGZsGPXbuP1Xy4CAF4c6oeXh3WutbJ0a2wJqsYgREREZKKu5RRiTuRZVKoFRPRyx7+Hd4ZEIsEAX2exS2s27BojIiIyQbcKyzBj42kUllair5cj1k7oCYmk9bb81IdBiIiIyMSUVqjwzOYzSLtTAi9nG3w5rQ+sLMzELksUDEJEREQmRK0WsGDbecTezIPc2gIbZvSDs51M7LJEwyBERERkQv6zPxG/XsyEhZkEX0ztA18XO7FLEhWDEBERkYn48XQq1v2VBABYM76nSQ2Krg+DEBERkQk4du02lvwSBwCYN9QPj/XpIHJFLQODEBERUSt3Nft/0+TH9HbH/OGdxS6pxWAQIiIiasVuFZZh5qb/TZN/9zHTnCZfHwYhIiKiVuqf0+S/eqqvyU6Trw+DEBERUSukVgv497ZYzTT5jTP6wcnWUuyyWhwGISIiolZo7b5E/HYxCxZmEnw5rQ98THyafH0YhIiIiFqZH6JT8cXh/02T7+/DafL1YRAiIiJqRY5evY0lO/9/mnxoJ06Tvw8GISIiolbiSnYh5kaehap6mvywTmKX1OIxCBEREbUCtwrLMHPjaRSWVaJfR9PdTV5XDEJERERGrqRchac3n0F6Xgk6Otvgy2l9ITPnNHltMAgREREZsepp8udv5qGNTdVu8pwmrz0GISIiIiP27r4E/B73/9Pkp3KavK7MxS6AiIiItKNSC4hOViKnsBSu9la4fvsuvjx8HQCwdkJPPMBp8jpjECIiIjICUXGZWLEnHpn5pbXOvRTaCeMCOU2+MRiEiIiIWriouEzMjTwHoZ7zXd3sm7We1oRjhIiIiFowlVrAij3x9YYgCYC39sZDpa7vCmoIgxAREVELFp2srLM7rJoAIDO/FNHJyuYrqhVhECIiImrBcgrrD0GNuY5qYhAiIiJqwVztrfR6HdXEIERERNSCBXs7wcVeVu95CQCF3ArB3k7NV1QrYjRBKCIiAp6enrCysoJCocC0adOQkZGh1b2CIGDkyJGQSCTYuXOnYQslIiLSI7UgwE5W93YZ1TuJLQ/3h5mU+4o1htEEoSFDhmDbtm1ITEzE9u3bkZSUhAkTJmh170cffcSN54iIyCj9Z18ikm8Xw8pCWqtlyE1uhXVTgxAWoBCpOuNnNOsIzZ8/X/NnLy8vLF68GGPHjkVFRQUsLCzqvS82Nhbvv/8+zpw5A4WCPyhERGQ8DiVk46u/q1aO/mhiIIb7t6uxsnSwtxNbgprIaILQvZRKJbZs2YKQkJAGQ1BxcTEmT56Mzz//HG5ublp97rKyMpSVlWk+LigoaHK9REREusrIK8G/t50HAMwI6YiwgKr3sQG+3EZDn4ymawwAFi1aBFtbWzg7OyM1NRW7du1q8Pr58+cjJCQEY8aM0fprrF69GnK5XPPy8PBoatlEREQ6qVCpMW9rDPKKK9CjvRyvjeoqdkmtlqhBaPHixZBIJA2+EhISNNcvXLgQMTEx2L9/P8zMzPDUU09BEOpeSXP37t04dOgQPvroI51qeu2115Cfn6953bx5synfIhERkc4++OMKzqTcgb3MHJ9NDoTMvO7B0tR0EqG+JNEMbt26hdzc3Aav8fHxgaWlZa3jaWlp8PDwwPHjxzFgwIBa519++WV88sknkEr/l/VUKhWkUikGDRqEv/76S6saCwoKIJfLkZ+fDwcHB63uISIiaqy/EnMwY+NpAMDnk4MwuifHtzaGtu/foo4RcnFxgYuLS6PuVavVAFBjPM+9Fi9ejKeffrrGsR49euDDDz9EeHh4o74mERGRIWXll2rGBU3t78kQ1AyMYrD0qVOncPr0aQwcOBCOjo5ISkrCsmXL4Ovrq2kNSk9PR2hoKDZv3ozg4GC4ubnVOUDa09MT3t7ezf0tEBERNajy/8cFKYvK4a9wwNLR/mKXZBKMYrC0jY0NduzYgdDQUHTp0gWzZ89Gz549cfjwYchkVWsqVFRUIDExEcXFxSJXS0REpLuPDlxF9A0lbC3N8PmUIFhZcFxQczCKFqEePXrg0KFDDV7TsWPHegdOVxNxOBQREVG9jly9hc//ugYAWP1YT3i3tRW5ItNhFC1CRERErVVOQSle/iEWggBMCvZERC93sUsyKQxCREREIlGpBbz0Qyxyi8rR1c0ey8M5Lqi5MQgRERGJ5JODV3Hiei5sOC5INAxCREREIjh+7TY+OXQVALBqXAB8XexErsg0MQgRERE1s1uFZXjpx6pxQU/07YBxgR3ELslkMQgRERE1I5VawPwfY3GrsAyd29lhRUSA2CWZNAYhIiKiZvT5n9dw9NptWFuY4fPJQbC25LggMTEIERERNZOT13Px0YErAICVYwPQqZ29yBURgxAREVEzuH23DPO2xkAtAI8FdcCEPhwX1BIwCBERERmY+v/HBeUUlsHP1Q4rx3YXuyT6fwxCREREBrbucBKOXL0NKwspPp8cBBtLo9jhyiQwCBERERnQ6RtKfPBH1bigFRHd0cWN44JaEgYhIiIiA1EWlePF72OgUgsY29sdT/T1ELsk+gcGISIiIgNQqwUs2BaLrIJS+LS1xdvjekAikYhdFv0DOymJiIj0QKUWEJ2sRE5hKVztrRB78w7+TLwFS3MpPpscBDsZ33JbIv5fISIiaqKouEys2BOPzPzSWufeDO8Of3cHEaoibTAIERERNUFUXCbmRp6DUM95RxuLZq2HdMMxQkRERI2kUgtYsSe+3hAkAfDW3nio1PVdQWJjECIiImqk6GRlnd1h1QQAmfmliE5WNl9RpBMGISIiokbKKaw/BDXmOmp+DEJERESN5GpvpdfrqPkxCBERETVSsLcTFPL6Q44EgEJuhWBvp+YrinTCIERERNRIZlIJnhnkU+e56qUTl4f7w0zKhRRbKgYhIiKiRqpUqbEzNh0AIDOv+ZbqJrfCuqlBCAtQiFEaaYnrCBERETXSl39fx4W0fDhYmSPq5cFIyS3WrCwd7O3EliAjwCBERETUCIlZhfj4wFUAwPLw7nBvYw33NtYiV0W6YtcYERGRjipVaiz8+TzKVWoM7eqK8UHtxS6JGolBiIiISEf3domtHs9d5Y0ZgxAREZEO/tkl1s6BawQZMwYhIiIiLbFLrPVhECIiItISu8RaHwYhIiIiLbBLrHViECIiIrqPe7vEQtkl1qowCBEREd3HvV1i77BLrFVhECIiImoAu8RaNwYhIiKierBLrPVjECIiIqoHu8RaPwYhIiKiOiRmFeKjA1cAsEusNWMQIiIi+ofqLrEKlcAusVZO593ny8rKcOrUKaSkpKC4uBguLi4IDAyEt7e3IeojIiJqduwSMx1aB6Fjx47h448/xp49e1BRUQG5XA5ra2solUqUlZXBx8cHzz77LObMmQN7e3tD1kxERGQw93aJvRnBLrHWTquusYiICEycOBEdO3bE/v37UVhYiNzcXKSlpaG4uBhXr17F0qVLcfDgQXTu3Bl//PGHoesmIiLSu0qVGq/89L8usXGB7BJr7bRqERo9ejS2b98OCwuLOs/7+PjAx8cH06dPR3x8PDIzM/VaJBERUXP48u/ruJjOLjFTIhEEQRC7iJasoKAAcrkc+fn5cHBwELscIiIykMSsQjz66RFUqAR88EQvjA/qIHZJ1ATavn/rPFj6Xnfv3oVara5xjGGBiIiMzb1dYsO6sUvMlOg8fT45ORmjR4+Gra0t5HI5HB0d4ejoiDZt2sDR0dEQNRIRERnUvV1iq8axS8yU6NwiNHXqVAiCgA0bNqBdu3b8YSEiIqPGWWKmTecgdP78eZw9exZdunQxRD1ERETNpoJdYiZP566xfv364ebNm4aohYiIqFl9de8sMXaJmSSdW4S++eYbzJkzB+np6QgICKg1pb5nz556K46IiMhQ/tkl5souMZOkcxC6desWkpKSMHPmTM0xiUQCQRAgkUigUqn0WiAREZG+sUuMquncNTZr1iwEBgbixIkTuH79OpKTk2v811AiIiLg6ekJKysrKBQKTJs2DRkZGQ3e8/DDD0MikdR4zZkzx2A1EhFRy6VSCziRlItdselY8stFdokRgEa0CKWkpGD37t3w8/MzRD31GjJkCF5//XUoFAqkp6fjlVdewYQJE3D8+PEG73vmmWfw1ltvaT62sbExdKlERNTCRMVlYsWeeGTml9Y4Pj6oPbvETJzOQWjo0KE4f/58sweh+fPna/7s5eWFxYsXY+zYsaioqKh36w+gKvi4ubk1R4lERNQCRcVlYm7kOdS1jcK3x1PQ38cZYQGKZq+LWgadg1B4eDjmz5+PixcvokePHrVCSEREhN6Kq49SqcSWLVsQEhLSYAgCgC1btiAyMhJubm4IDw/HsmXL2CpERGQiVGoBK/bE1xmCqq3YE4/h/m4wk7J7zBTpHISqx9jc291UzdCDpRctWoTPPvsMxcXF6N+/P/bu3dvg9ZMnT4aXlxfc3d1x4cIFLFq0CImJidixY0e995SVlaGsrEzzcUFBgd7qJyKi5hWdrKzVHXYvAUBmfimik5UY4OvcfIVRi6HzYGm1Wl3vS9cQtHjx4lqDmf/5SkhI0Fy/cOFCxMTEYP/+/TAzM8NTTz2FhvaMffbZZzFixAj06NEDU6ZMwebNm/HLL78gKSmp3ntWr14NuVyueXl4eOj0PRERUcuRU1h/CGrMddT6iLr7/K1bt5Cbm9vgNT4+PrC0tKx1PC0tDR4eHjh+/DgGDBig1dcrKiqCnZ0doqKiMGLEiDqvqatFyMPDg7vPExEZoRNJuZj09cn7Xrf1mf5sEWplDLr7/OHDh/Hee+/h8uXLAAB/f38sXLgQgwYN0unzuLi4wMXFpTElaHa9vze03E9sbCwAQKGof1CcTCaDTCZrVE1ERNSyBHs7wdVehpzCut8rJADc5FYI9nZq3sKoxdC5aywyMhLDhg2DjY0N5s2bh3nz5sHa2hqhoaH4/vvvDVEjTp06hc8++wyxsbFISUnBoUOHMGnSJPj6+mpag9LT09G1a1dER0cDAJKSkrBy5UqcPXsWN27cwO7du/HUU09h8ODBXP2aiMhEmEkl8HO1q/Nc9dDo5eH+HChtwnRuEVq1ahXWrl1bYzr7vHnz8MEHH2DlypWYPHmyXgsEqqbA79ixA8uXL0dRUREUCgXCwsKwdOlSTetNRUUFEhMTUVxcDACwtLTEgQMH8NFHH6GoqAgeHh547LHHsHTpUr3XR0RELdPpG0ocT6oaguFsa4nconLNOTe5FZaH+3PqvInTeYyQTCbDpUuXaq0jdO3aNQQEBKC0tHUNONO2j5GIiFqWCpUaj35yFInZhXiynwdWjeuB6GQlcgpL4Wpf1R3GlqDWy2BjhDw8PHDw4MFaQejAgQOcYUVERC3GhqPJSMwuhJOtJRaFdYWZVMIB0VSLzkFowYIFmDdvHmJjYxESEgIAOHbsGDZt2oSPP/5Y7wUSERHpKu1OMT46cBUA8NrIrnC0rT37mAhoRBCaO3cu3Nzc8P7772Pbtm0AgG7duuHHH3/EmDFj9F4gERGRrt7cHY+SChWCvZ0woU8HscuhFqxR0+fHjRuHcePG6bsWIiKiJtt/KQsHLmfDXCrB22MDuLM8NUjn6fM+Pj51LoKYl5cHHx8fvRRFRETUGMXllVixJx4A8MxgH3RuZy9yRdTS6RyEbty4UedWGmVlZUhPT9dLUURERI3x8cGrSM8rQQdHa8wb2knscsgIaN01tnv3bs2f9+3bB7lcrvlYpVLh4MGD6Nixo16LIyIi0lZiViHWH0kGAKyI6A5rSzORKyJjoHUQGjt2LICqHeanT59e45yFhQU6duyI999/X6/FERERaUOtFrB050VUqgWM6N4Ood3aiV0SGQmtg1D13l7e3t44ffo02rZta7CiiIiIdPHz2TScvnEHNpZmWB7eXexyyIjoPGssOTnZEHUQERE1irKoHO/8XrUJ+PxhneHexlrkisiYNGr6fFFREQ4fPozU1FSUl5fXODdv3jy9FEZERKSNNb9fRl5xBbq62WPGgx3FLoeMjM5BKCYmBqNGjUJxcTGKiorg5OSE27dvw8bGBq6urgxCRETUbE7fUGLbmTQAwKpxAbAw03kyNJk4nX9i5s+fj/DwcNy5cwfW1tY4efIkUlJS0KdPH7z33nuGqJGIiKiWCpUaS365CACYFOyBPl5OIldExkjnIBQbG4sFCxZAKpXCzMwMZWVl8PDwwNq1a/H6668bokYiIqJa1h9NxpXsu5pNVYkaQ+cgZGFhAam06jZXV1ekpqYCAORyOW7evKnf6loplVrAiaRc7IpNx4mkXKjUgtglEREZlbQ7xfj4/zdVfX1UN7Sx4aaq1Dg6jxEKDAzE6dOn0alTJzz00EN44403cPv2bXz33XcICAgwRI2tSlRcJlbsiUdmfqnmmEJuheXh/ggLUIhYGRGR8bh3U9XHgtqLXQ4ZMZ1bhN555x0oFFVv2KtWrYKjoyPmzp2LW7du4auvvtJ7ga1JVFwm5kaeqxGCACArvxRzI88hKi5TpMqIiIzHvZuqruKmqtREWrcIFRcXw8bGBn379tUcc3V1RVRUlEEKa21UagEr9sSjrk4wAYAEwIo98Rju7wYzKf9SExHVpaisEm/uvgSgalPVTtxUlZpI6xahtm3b4tFHH8VXX32FrKwsQ9bUKkUnK2u1BN1LAJCZX4roZGXzFUVEZGQ+OXgVGfml3FSV9EbrIJSQkIARI0Zg27Zt6NixIx544AGsWrUKFy9eNGR9rUZOYf0hqDHXERGZmoSsAqw/WrW7wVtjuKkq6YfWQcjT0xMvvvgiDhw4gOzsbLz88su4ePEiBg0aBB8fH7z88ss4dOgQVCqVIes1Wq72Vnq9jojIlKjVApb+EqfZVHVoV26qSvrRqCU45XI5Jk2ahB9++AG3bt3Cl19+CZVKhZkzZ8LFxQVbtmzRd51GL9jbCQq5Feob/SNB1eyxYG8uCEZE9E8/n03DmRRuqkr6p3UQiouLq/O4hYUFhg8fjk8//RQpKSk4ePAgOnfurLcCWwszqQTLw/0BoM4wJABYHu7PgdJERP/ATVXJkLQOQj179sQDDzyAr7/+GoWFhfVeFxgYiH79+umluNYmLECBdVOD4Cav3f1lZS5F345sDSIi+qfVv3FTVTIcrYPQ4cOH0b17dyxYsAAKhQLTp0/HkSNHDFlbqxQWoMDRRUOx9Zn++PjJ3vj+6QcQ4O6A0kq1ZpVUIiKqEp2sxE9nuakqGY7WP1GDBg3Chg0bkJmZiU8//RQ3btzAQw89hM6dO+Pdd9/llHodmEklGODrjDG92yPEry2WjK7qMvs+OhXXcu6KXB0RUctQoVJj6U5uqkqGpXO0trW1xcyZM3H48GFcuXIFjz/+OD7//HN4enoiIiLCEDW2egN8nTGsmytUagHvRiWIXQ4RUYvATVWpOTSpjdHPzw+vv/46li5dCnt7e/z666/6qsvkLB7ZFWZSCf6Iz8bJ67lil0NEJCpuqkrNpdFB6O+//8aMGTPg5uaGhQsXYvz48Th27Jg+azMpfq72eLKfBwDgnd8uQ80d6YnIxKjUAk4k5WJXbDrmbY1FSYUKD3BTVTIwnXafz8jIwKZNm7Bp0yZcu3YNISEh+OSTT/DEE0/A1tbWUDWajJeHdcbOmHRcSMvHngsZGNObf/mJyDRExWVixZ74WlsRjejuxk1VyaC0bhEaOXIkvLy88Omnn2LcuHG4fPkyjh49ipkzZzIE6YmLvQxzH/YFAKyNSkRpBVfpJqLWLyouE3Mjz9W5H+PKvfGIissUoSoyFVoHIQsLC/z8889IS0vDu+++iy5duhiyLpM1e6AP3ByskJ5Xgs0nbohdDhGRQanUAlbsiUdDgwFW7ImHisMFyEC0DkK7d+/GmDFjYGZWtcndtWvXsG/fPpSUlAAABIE/pPpgbWmGBY9Urcz96aFruFNULnJFRESGE52srLMlqJoAIDO/FNHJyuYrikyKzoOlc3NzERoais6dO2PUqFHIzKxqspw9ezYWLFig9wJN0figDuimcEBhaSU+PXRN7HKIiAwmp7D+ENSY64h0pXMQmj9/PiwsLJCamgobGxvN8YkTJyIqKkqvxZkqM6kEr4+qWjPju5M3cON2kcgVEREZhqt97S2HmnIdka50DkL79+/Hu+++iw4dOtQ43qlTJ6SkpOitMFM3qJMLHursggqVgLX7uMgiEbVOwd5OaOcgq/e8BIBCboVgb64qTYahcxAqKiqq0RJUTalUQiar/4eZdPfaqK6QSoDfLmbhbModscshItI7M6kEvT3a1HmuetL88nB/mEk5hZ4MQ+cgNGjQIGzevFnzsUQigVqtxtq1azFkyBC9Fmfquro54PE+VYssrvo1ngPSiajVuXG7CH8m3AIAONpY1DjnJrfCuqlBCAtQiFEamQidFlQEgLVr1yI0NBRnzpxBeXk5Xn31VVy6dAlKpZIrSxvAvx/pjN3nM3AuNQ9RcVkY2YO/EIio9XhrbzzKVWoM7uyCDdP74vSNO8gpLIWrfVV3GFuCyNB0bhEKCAjAlStXMHDgQIwZMwZFRUUYP348YmJi4Ovra4gaTVo7Bys8M9gHALAmKgHllWqRKyIi0o9DCdk4lJADCzMJlof7w9xMigG+zhjTuz0G+DozBFGz0LlFCADkcjmWLFmi71qoHs8N9sH3p1KRkluMLadSMPNBb7FLIiJqktIKFVbsiQcAzHrQG74udiJXRKZKqxah1NRUnT5penp6o4qhutnKzDWLLH588CrySypEroiIqGnWH01GSm4xXO1leDG0k9jlkAnTKgj169cPzz33HE6fPl3vNfn5+fj6668REBCA7du3661AqvJ4nw7o5GqHvOIK/PdPLrJIRMYrI68En/3/YrGvj+oGO1mjOieI9EKrn774+HisWrUKw4cPh5WVFfr06QN3d3dYWVnhzp07iI+Px6VLlxAUFIS1a9di1KhRhq7b5JibSfH6qG6Yuek0Nh6/gan9veDhVHsZAyKilu6d3y6jpEKFfh0dMaa3u9jlkInTqkXI2dkZH3zwATIzM/HZZ5+hU6dOuH37Nq5evQoAmDJlCs6ePYsTJ04wBBnQw11cEOLrjPJKNd7bnyh2OUREOjuRlIu9FzIhlQBvRnSHRMIB0SQuicDFaRpUUFAAuVyO/Px8ODg4iF0O4tLzEf7ZUQgCsPuFB9GzQxuxSyIi0kqlSo3RnxxFYnYhpvb3xNtje4hdErVi2r5/6zx9nsQV0F6OcYHtAQBv743HiaTb2BWbjhNJuVCpmWmJqOX67mQKErML4WhjgVce6SJ2OUQAGjl9nsT1yiNdsDs2A9E37mDS16c0xxVyKywP9+cqrETU4ty+W4YP/rgCAHhlRBe0sbEUuSKiKmwRMkIX0vJQWUfrT1Z+KeZGnkNUXKYIVRER1e8/UYkoLK1EQHsHPNnPU+xyiDQYhIyMSi1oFiH7p+potGJPPLvJiKjFiL2Zh21nbwIAVkR054rR1KI0avd5Ek90shKZ+aX1nhcAZOaXIjpZ2XxFERHVQ60WsHxXHAQBGB/YHn28nMQuiagGnYNQu3btMGvWLBw9etQQ9dB95BTWH4Iacx0RkSH9fDYN59PyYSczx+KRXcUuh6gWnYNQZGQklEolhg4dis6dO2PNmjXIyMgwRG01REREwNPTE1ZWVlAoFJg2bZpWX/fEiRMYOnQobG1t4eDggMGDB6OkpMTg9RqKq72VXq8jIjKU/JIKvBuVAAB4KbQTXB34e4laHp2D0NixY7Fz506kp6djzpw5+P777+Hl5YVHH30UO3bsQGVlpSHqxJAhQ7Bt2zYkJiZi+/btSEpKwoQJExq858SJEwgLC8MjjzyC6OhonD59Gi+88AKkUuMdGhXs7QSF3Ar19bBLUDV7LNibzc9EJK6PDlxBblE5fF1sMT2ko9jlENVJLwsqfvrpp1i4cCHKy8vRtm1bzJkzB4sXL4aNjeG2gNi9ezfGjh2LsrIyWFhY1HlN//79MXz4cKxcubLRX6elLagIAFFxmZgbeQ7A/wZI3+uLqUGcQk9EokrMKsSoT45ApRbw3exgDOrkInZJZGIMvqBidnY21q5dC39/fyxevBgTJkzAwYMH8f7772PHjh0YO3ZsYz/1fSmVSmzZsgUhISH1hqCcnBycOnUKrq6uCAkJQbt27fDQQw/dd2xTWVkZCgoKarxamrAABdZNDYKbvHYzs72VOUL82opQFRFRFUEQsHx3HFRqASO6t2MIohZN5wUVd+zYgY0bN2Lfvn3w9/fH888/j6lTp6JNmzaaa0JCQtCtWzd91gkAWLRoET777DMUFxejf//+2Lt3b73XXr9+HQDw5ptv4r333kPv3r2xefNmhIaGIi4uDp06darzvtWrV2PFihV6r13fwgIUGO7vhuhkJXIKS+FkY4k3dsUhObcYHx+4imWP+otdIhGZqF8vZuLkdSVk5lIsHc3fRdSy6dwiNHPmTLi7u+PYsWOIjY3FCy+8UCMEAYC7uzuWLFly38+1ePFiSCSSBl8JCQma6xcuXIiYmBjs378fZmZmeOqpp1Bfz55arQYAPPfcc5g5cyYCAwPx4YcfokuXLtiwYUO9Nb322mvIz8/XvG7evKnFUxGHmVSCAb7OGNO7PQZ1dsGKMQEAgE3HbyAxq1Dk6ojIFBWXV2LVr5cBAHMf9oWHk+GGSBDpg84tQpmZmfcd+2NtbY3ly5ff93MtWLAAM2bMaPAaHx8fzZ/btm2Ltm3bonPnzujWrRs8PDxw8uRJDBgwoNZ9CkXVGBl//5r/GunWrRtSU1Pr/XoymQwymey+tbdEgzu7IKy7G6IuZWH57jhsfaY/d3Ymomb13z+TkJlfig6O1pjzkK/Y5RDdl85BqLKyss5xMxKJBDKZDJaW2u8f4+LiAheXxvUdV7f4lJWV1Xm+Y8eOcHd3R2JiYo3jV65cwciRIxv1NY3B0ke74a8rOTh5XYk9FzIR0ctd7JKIyETcuF2Er/6uGpawdLQ/rCzMRK6I6P507hpr06YNHB0da73atGkDa2treHl5Yfny5Zqgog+nTp3CZ599htjYWKSkpODQoUOYNGkSfH19Na1B6enp6Nq1K6KjowFUBbOFCxfik08+wc8//4xr165h2bJlSEhIwOzZs/VWW0vTwdEG/3rYDwCw6td4FJUZZjkDIqJ/Wrk3HuUqNQZ1aosR3duJXQ6RVnRuEdq0aROWLFmCGTNmIDg4GAAQHR2Nb7/9FkuXLsWtW7fw3nvvQSaT4fXXX9dLkTY2NtixYweWL1+OoqIiKBQKhIWFYenSpZpurIqKCiQmJqK4uFhz38svv4zS0lLMnz8fSqUSvXr1wh9//AFf39bdXPvMYB/8fC4NKbnF+OTQVbw2Uv8D14mI7vVnQg4OJuTAXCrB8vDu7JYno6HzOkKhoaF47rnn8MQTT9Q4vm3bNnz55Zc4ePAgvvvuO6xatarGQGdj1RLXEdLGoYRszNp0BhZmEvz+0mD4udqJXRIRtVJllSqM+PBv3MgtxrODffD6KP7ji8RnsHWEjh8/jsDAwFrHAwMDceLECQDAwIEDGxyQTIY3tGs7hHZ1RYVKwIo9l+qdXUdE1FTrjybjRm4xXOxleHGon9jlEOlE5yDk4eGB9evX1zq+fv16eHh4AAByc3Ph6OjY9OqoSd4I94eluRRHrt7GvktZYpdDRK1QZn4JPj14DQDw2siusLeqe5FbopZK5zFC7733Hh5//HH8/vvv6NevHwDgzJkzSEhIwM8//wwAOH36NCZOnKjfSklnXs62mDPYB58cuoaVey/joc6usLbkLA4i0p93fktASYUKfb0cMS6wvdjlEOmsUXuN3bhxA19++aVmanqXLl3w3HPPoWPHjvquT3TGOkaoWkm5CsM+OIz0vBK8ONQPCx7pInZJRNRKnLyeiye/OgmJBNjzwkAEtJeLXRKRhrbv3zq1CFVUVCAsLAxffPEFVq9e3eQiyfCsLc2w7FF/zIk8iy8PX8djQR3Qsa2t2GURkZGrVKnx5u5LAIDJwZ4MQWS0dBojZGFhgQsXLhiqFjKQEd3bYXBnF5Sr1Hhrb7zY5RCRkVKpBZxIysWu2HS8/etlJGQVoo2NBV5hSzMZMZ0HS0+dOrXOwdLUckkkErwZ7g8LMwkOJeTgQHy22CURkZGJisvEwHcPYdLXJ/HSD7HYdPwGACAswA2OttrvKEDU0jRqi40NGzbgwIED6NOnD2xta3azfPDBB3orjvTHx8UOTw/ywbq/krBi7yUM7NSWy98TkVai4jIxN/Ic6hpQ+mP0TTzc2QVhAYpmr4tIH3QOQnFxcQgKCgJQtW/XvbiSaMv2whA/7IxJx01lCb48fB0vDeskdklE1MKp1AJW7ImvMwRVW7EnHsP93WAm5XsAGR+dg9Cff/5piDqoGdjKzLFkdDe88H0M/vvXNYwPag8PJxuxyyKiFiw6WYnM/NJ6zwsAMvNLEZ2sxABf5+YrjEhPdB4jVO3atWvYt28fSkpKAIArFxuJ0T0UCPF1RlmlGis5cJqI7iOnsP4Q1JjriFoanYNQbm4uQkND0blzZ4waNQqZmZkAgNmzZ2PBggV6L5D0SyKRYEVEd5hLJdgfn42/EnPELomIWjBXeyu9XkfU0ugchObPnw8LCwukpqbCxuZ/3SoTJ05EVFSUXosjw+jUzh4zQjoCqOrbLy6v1EyJPZGUC5WarXtEVCXY2wkKef0hRwJAIbdCsLdT8xVFpEc6jxHav38/9u3bhw4dOtQ43qlTJ6SkpOitMDKsl4Z1wq7zGUi+XYQHVh1EYVml5pxCboXl4f6cBUJEMJNKMDnYE+//caXWueqh0cvD/TlQmoyWzi1CRUVFNVqCqimVSshkMr0URYZnb2WBR3u6AUCNEAQAWfmlmBt5DlFxmWKURkQtSIVKjb0Xqn4X/HOvQje5FdZNDeI/msio6dwiNGjQIGzevBkrV64EUDXmRK1WY+3atRgyZIjeCyTDUKkF/B5X98KKAqr+pccpsUS04WgyErML4WRriT/mD8aV7LvIKSyFq31Vdxh/P5Cx0zkIrV27FqGhoThz5gzKy8vx6quv4tKlS1AqlTh27JghaiQDiE5WIotTYomoARl5JfjowFUAwGsju8LZToYBdmz5p9ZF566xgIAAXLlyBQMHDsSYMWNQVFSE8ePHIyYmBr6+voaokQyAU2KJ6H5W7LmEkgoV+nV0xGNBHe5/A5ER0rlFCADkcjmWLFmi71qoGXFKLBE15FBCNvZdyoaZVIK3x/aAlF1g1Eo1Kgjl5eUhOjoaOTk5UKvVNc499dRTeimMDKt6SmxWfmmdS+dLUDUQklNiiUxPSbkKb+y6BACYPdAbXdzsRa6IyHB0DkJ79uzBlClTcPfuXTg4ONTYX0wikTAIGQkzqQTLw/0xN/IcJECdYYhTYolM0+d/XkPanRK4y63wUij3JKTWTecxQgsWLMCsWbNw9+5d5OXl4c6dO5qXUqk0RI1kIGEBCqybGgS3OhZLG91TwSmxRCYo6dZdfPl3EgDgjfDusJU1quOAyGjo/BOenp6OefPm1bmWEBmfsAAFhvu7ITpZiZzCUlzLuYtPD13DvktZuJJdiM7t2CROZCoEQcCynXGoUAkY0sUFI7q3E7skIoPTuUVoxIgROHPmjCFqIZGYSSUY4OuMMb3b49/DO2NYt3aoUAl49ecL3G6DyITsPp+B40m5kJlLsSIioMbQB6LWSucWodGjR2PhwoWIj49Hjx49YGFhUeN8RESE3oqj5ieRSPD22ACcup6L2Jt52HgsGU8P8hG7LCIysILSCqzcexkA8OJQP3g6s9WfTINEEASd/skvldbfiCSRSKBSqZpcVEtSUFAAuVyO/Px8ODg4iF1Os9kanYrXdlyElYUU+19+iL8UiVq55bvi8O2JFPi42OL3lwZBZm52/5uIWjBt37917hpTq9X1vlpbCDJlT/bzwAAfZ5RWqLF4xwXomJeJyIhcTMvHdyerNs1+e0wAQxCZFJ2DEJkGiUSCNY/1gJWFFMeTcvHj6Ztil0REBqBSC1i68yLUAjCmtztC/NqKXRJRs9I6CI0aNQr5+fmaj9esWYO8vDzNx7m5ufD399drcSQuL2dbvPJIFwDAql8vN7g3GREZp++jU3E+LR/2MnMsGd1N7HKImp3WQWjfvn0oKyvTfPzOO+/UWDeosrISiYmJ+q2ORDfzQW/08miDwrJKLN0Zxy4yolbkVmEZ1kYlAABeGdGFW+qQSdI6CP3zDZBviKbBTCrB2sd6wsJMggOXs7H3QqbYJRGRnqz+7TIKSysR0N4BU/t7iV0OkSg4Rojuq4ubPf41xA8A8ObuS1AWlYtcERE11YmkXOyISYdEAqwa24Pb6ZDJ0joISSSSWotrcbEt0/H8w37o0s4euUXleGvPJbHLIaImKK9UY9muOADAlAc80cujjbgFEYlI6wUVBUHAjBkzIJPJAAClpaWYM2cObG1tAaDG+CFqfSzNpXh3Qk+M/+8x7IzNQERvdwztyuX3iYzRN0ev41rOXbS1s8TCR7qKXQ6RqLRuEZo+fTpcXV0hl8shl8sxdepUuLu7az52dXXlzvOtXG+PNpg90BsAsOSXOBSWVohcERHp6qayGJ8cvAoAeH1UN8htLO5zB1HrpnWL0MaNGw1ZBxmJfw/vgv3x2UjJLcaa3xOwalwPsUsiIh2s2BOP0go1+vs4YVxge7HLIRIdB0uTTqwtzbB6fFX42XIqFSev54pcERFp64/4bBy4nA1zadWeghznScQgRI0Q4tsWk4I9AQCLt19ASTm3ViFq6YrLK/Hm7qqJDs8M9oGfq73IFRG1DAxC1CivjeoKNwcr3MgtxkcHrohdDhHdx6eHriE9rwTt21hj3tBOYpdD1GIwCFGjOFhZ4O2xAQCAr49cx4W0PHELIqJ6Xc0uxNd/XwcArIjoDmtLbqpKVI1BiBptmH87RPRyh1oAXv25qovsRFIudsWm40RSLlRqrj5OJDZBELB0Zxwq1QKGdWuHYf5c9oLoXlrPGiOqy/Jwfxy9dhsJWYXot+oA7pZVas4p5FZYHu6PsACFiBUSmbZfYtJxKlkJawszvBnBjbGJ/oktQtQkznYyjO3tDgA1QhAAZOWXYm7kOUTFcX8yIjHkF1dg1a+XAQDzQjuhg6ONyBURtTxsEaImUakF/BaXVec5AYAEVeuWDPd3415GRM1ApRYQnaxETmEp9pzPQG5ROTq52mkWQyWimhiEqEmik5XIyi+t97wAIDO/FNHJSgzwdW6+wohMUFRcJlbsiUfmP/5Ohvdyh6U5OwCI6sK/GdQkOYX1h6DGXEdEjRMVl4m5kedqhSAA+PCPK+yiJqoHgxA1iau9lV6vIyLdqdQCVuyJR0PzNFfsiedMTqI6MAhRkwR7O0Eht0J9o38kqJo9Fuzt1JxlEZmU6GRlnS1B1e7toiaimhiEqEnMpBIsD6+akltXGBJQNcWeA6WJDIdd1ESNxyBETRYWoMC6qUFwk9fu/rK1NEOgp6MIVRGZDnZREzUeZ42RXoQFKDDc300zbdfRxgKrf0vA5axC/HtbLL6b9QCkbBUiMohgbye0tbPE7bvldZ6XAHBjFzVRnYymRSgiIgKenp6wsrKCQqHAtGnTkJGRUe/1N27cgEQiqfP1008/NWPlpsNMKsEAX2eM6d0egzu74tPJQbCykOLYtVx8c/S62OURtVqCIMBWVve/a6v/+cEuaqK6GU0QGjJkCLZt24bExERs374dSUlJmDBhQr3Xe3h4IDMzs8ZrxYoVsLOzw8iRI5uxctPl52qH5eHdAQD/2ZeIi2n5IldE1Dr9968kpOQWw9pCCld7WY1zbnIrrJsaxK1uiOohEQTBKOdT7t69G2PHjkVZWRksLCy0uicwMBBBQUFYv3691l+noKAAcrkc+fn5cHBwaGy5JksQhKptNi5lwaetLfbOGwgbS/bIEunLxbR8jPvvMVSqBXz8ZG882tNd00Xtal/VHcaWIDJF2r5/G+U7klKpxJYtWxASEqJ1CDp79ixiY2Px+eefN3hdWVkZysrKNB8XFBQ0qVZTJ5FIsOaxHoi9mYfrt4vw1p54rHmsp9hlEbUKpRUqzN8Wi0q1gNE9FYjo5Q6JRMJV3Il0YDRdYwCwaNEi2NrawtnZGampqdi1a5fW965fvx7dunVDSEhIg9etXr0acrlc8/Lw8Ghq2SavjY0lPpjYCxIJ8MPpm/jtIle4JdKH/+xLxLWcu3Cxl+HtMQGQSNjyQ6QrUYPQ4sWL6x3QXP1KSEjQXL9w4ULExMRg//79MDMzw1NPPQVtevZKSkrw/fffY/bs2fe99rXXXkN+fr7mdfPmzSZ9j1QlxLct5j7kCwBYvP0CMvJKRK6IyLgdT7qN9UeTAQBrJ/SEo62lyBURGSdRxwjdunULubm5DV7j4+MDS8vaf8HT0tLg4eGB48ePY8CAAQ1+ju+++w6zZ89Geno6XFxcdKqRY4T0p0KlxoR1x3E+LR/B3k7Y+kx/jl0gaoSC0gqM/OgI0vNKMPkBT7wzrofYJRG1OEYxRsjFxUXnYFJNrVYDQI3xPPVZv349IiIiGv21SD8szKT4+MlAjP7kCKKTlfjicBL+NcRP7LKIjM5be+KRnlcCTycbLBnVTexyiIyaUYwROnXqFD777DPExsYiJSUFhw4dwqRJk+Dr66tpDUpPT0fXrl0RHR1d495r167h77//xtNPPy1G6fQPHdvaYsWYAADAB39cQUzqHZErIjIu+y5l4eezaZBIgA+e6FXv+kFEpB2jCEI2NjbYsWMHQkND0aVLF8yePRs9e/bE4cOHIZNVrZlRUVGBxMREFBcX17h3w4YN6NChAx555BExSqc6PBbUHuG93KFSC3jph1gUllaIXRKRUbh9twyv77gIAHhusC/6duRK0URNZbTrCDUXjhEyjPySCoz6uGqMw/jA9vhgYm+xSyJq0QRBwLPfncUf8dno6maPXS88CJm5mdhlEbVY2r5/G0WLELU+cmsLfPxkb0glwI6YdOyKTRe7JKIW7eezafgjPhsWZhJ8OLE3QxCRnjAIkWj6dnTCi0M7AQCW/hKHm8ri+9xBZJrS7hRjxZ54AMC/h3dBNwVbp4n0hUGIRPXiUD/08XJEYVklXvohBpUqtdglEbUoarWAV346j7tllejj5YhnB/uIXRJRq8IgRKIyN5Pio4m9YS8zx7nUPHxy6JrYJRG1KBuOJePkdSVsLM3wwRO9uPYWkZ4xCJHoPJxs8Pa4qin1nx26ipPXc3EiKRe7YtNxIikXKjXH85NpuppdiLX7EgEAS0f7w8vZVuSKiFofLkBBLcKY3u1x+Mot7DiXjslfn8S92Ucht8LycH+EBSjEK5ComVWo1Ji/LRbllWo83MUFk4K57yGRIbBFiFqMQZ2qVv7+ZwNQVn4p5kaeQ1QcN2sl0/HpoWuISy9AGxsLrH2sJzdUJTIQBiFqEVRqAWujEuo8V52LVuyJZzcZmYTYm3n4/M+q8XJvjw2Aq4OVyBURtV4MQtQiRCcrkZlfWu95AUBmfimik5XNVxSRCErKVfj3j7FQqQVE9HLHoz3dxS6JqFVjEKIWIaew/hDUmOuIjNW7UQm4frsI7RxkWPn/+/IRkeEwCFGL4GqvXdO/ttcRGaOjV29j0/EbAID/TOgFuY2FuAURmQAGIWoRgr2doJBboaHhoAq5FYK9uckktU75JRVY+PN5AMC0/l4Y3NlF5IqITAODELUIZlIJlof7A0C9YWi4fzsuJket1pu7LyEzvxQdnW3w2qiuYpdDZDIYhKjFCAtQYN3UILjJa3Z/2cmqlrv6IfomzqbcEaM0IoP67WImfolJh1QCfDCxN2wsucQbUXPh3zZqUcICFBju74boZCVyCkvham+Fvl6OeGHrOey7lI05kWex54WBtcISkTFRqQXNz7iluRSv77gAAHj+YT8EeTqKXB2RaZEIgsCFWRpQUFAAuVyO/Px8ODhwx2exFJVVYvx/jyMxuxA9O8ix7bkBsLIwE7ssIp1FxWVixZ74WstFdGhjjUOvPAxLczbUE+mDtu/f/BtHRsFWZo5vpveFo40FLqTlY/H2C2CGJ2MTFZeJuZHn6lwzKy2vBIcSskWoisi0MQiR0fBwssHnU4JgJpVgZ2wGvvr7utglEWlNpRawYk886ovvEnD1dCIxMAiRUQnxbauZXbYmKgF/JuaIXBGRdrh6OlHLxCBERmdafy9MCvaAIADzvo/BtZy7YpdEdF9cPZ2oZWIQIqMjkUiwIiIA/To6orCsEs9uPoP8kgqxyyJqEFdPJ2qZGITIKFmaS7Fuah+4y61w/XYR5m2N4dgKatGCvZ3Qxrr+LTMk4OrpRGJgECKj1dZOhq+e6gsrCykOX7mFd6MSxC6JqF6xN+/gblllneeq10tfHu7P1dOJmhmDEBm1gPZyvPd4LwDAV39fx45zaSJXRFTbTWUxnt18FpVqAb06yOHmULP7y01uhXVTgxAWoBCpQiLTxZWlyeg92tMdCZmF+OzPa1i84yJ8XOzQ26ON2GURAQAKSiswa9Np5BaVo7u7A7Y+2x8yc7Maq6cHezuxJYhIJGwRolbh38M7Y1i3diivVOO5784gp4Azb0h8lSo1Xvg+Bldz7qKdgwzrp/eDjaU5zKQSDPB1xpje7THA15khiEhEDELUKkilEnw4sRc6udohu6AMz353FqUVKrHLIhO3cm88/r5yC9YWZlg/vR/3yCNqgRiEqNWwt7LAN9P7Qm5tgdibeXj9l4uoVKlxIikXu2LTcSIplzPLqNl8e/wGvj2RAokE+HBibwS0l4tdEhHVgZuu3gc3XTU+R6/exvSN0VCpBdhbmaOw9H8zdRRyKywP9+egVDKovxJzMGvTaagFYPHIrpjzkK/YJRGZHG66SiZrYKe2GB/YHgBqhCAAyMovxdzIc4iKyxSjNDIBiVmFeOH7GKgF4PE+HfDcYB+xSyKiBjAIUaujUgs4cu12neeqmz+5uSUZwu27ZZi16TTullXiAW8nrBrXAxIJB0ITtWQMQtTqRCcrkcXNLamZlVao8OzmM0jPK0FHZxt8MbUPLM35K5aopePfUmp1uLklNTdBEPDqzxdwLjUPDlbmWD+jHxxtLcUui4i0wCBErQ43t6Tm9vHBq9h9PgPmUgm+mNoHvi52YpdERFpiEKJWJ9jbCQq5FRoamdHOQcbNLUkvdsWm46MDVwEAb48NQIhfW5ErIiJdMAhRq2MmlWB5uD8A1BuGrC3MUFxe9waYRNo6m3IHC3++AAB4drAPngz2FLkiItIVgxC1SmEBCqybGlRrJd+2djLYWJrhRm4xZmw8Xe9u4ET3c1NZjOe+O4PySjWGdWuHRWFdxS6JiBqBCyreBxdUNG4qtVBrc8vLmQWY/PVJFJRWol9HR2yaGQxbGfcfJu0VllZgwroTSMwuhL/CAT/NGcCfIaIWhgsqEgF1bm4Z0F6OyKcfgL2VOU7fuINZm06zm4y0VqlS48WtMUjMLoSrvQzrZ/RlCCIyYgxCZJJ6dmiDzbOCYSczx6lkJZ7+9gw3aSWtvP3rZfyVeAtWFlJ8M70vFHJrsUsioiZgECKTFejpiG9n9YOtpRmOJ+Ximc0MQ1STSi3U2LR30/FkbDp+AwDw0cTe6Nmhjaj1EVHTcYzQfXCMUOsXnazE9A3RKKlQ4eEuLvhyWh/IzM3ELotEFhWXiRV74pFZxyrlr4Z1wfMP+4lQFRFpi2OEiLQU7O2EDTP6wcpCir8Sb+H5yHMor1SLXRaJKCouE3Mjz9UZggDA29m2mSsiIkNhECICMMDXGeun94PMXIqDCTl44ftzqFAxDJkilVrAij3xqK+pXALgrb3ctJeotWAQIvp/D/q1xddP9YWluRT747Mxb2sMw5AJik5W1tsSBHDTXqLWhkGI6B6DO1eNEbI0k+L3uCzM/zEWlQxDJoWb9hKZFgYhon8Y0sUV66YGwcJMgr0XMrHgp/PsBjEh3LSXyLQwCBHVIbRbO3w+OQjmUgl2xWZg4c/nUV6prjGVmuGodcq9W9bgeQkAhdyKm/YStRKcPn8fnD5v2n6/mIkXtsZApRZgbWGGknvWGVLIrbA83B9hAQoRKyR9EQQBX/19Hat/T9AckwA1Bk1Xb+K7bmoQ/78TtXCcPk+kByN7KDDzwY4AUCMEAUBWfinmRp5DVFymCJWRPlWq1Fi6M04TgmaEdMR/J9fetNdNbsUQRNTKcIMcogao1AJ+vVB30BFQ1UKwYk88hvu7wUwqqfM6atnullXihe/P4a/EW5BIgGWj/TFroDcAYESAW61Ne/n/mah1YRAiaoAuU6kH+Do3X2GkF1n5pZi56TQuZxbAykKKj58MxIjubprz1Zv2ElHrZTRdYxEREfD09ISVlRUUCgWmTZuGjIyMBu/JysrCtGnT4ObmBltbWwQFBWH79u3NVDG1BpxK3XrFZxRg7OfHcDmzAG3tLPHjswNqhCAiMg1GE4SGDBmCbdu2ITExEdu3b0dSUhImTJjQ4D1PPfUUEhMTsXv3bly8eBHjx4/HE088gZiYmGaqmowdp1K3Tn8m5uDxL44jq6AUfq52+OX5B9HLo43YZRGRCIx21tju3bsxduxYlJWVwcLCos5r7OzssG7dOkybNk1zzNnZGe+++y6efvpprb4OZ42ZNpVawMB3DyErv7TeLRcA4Im+HfDWmABYWXCz1pZuy6kUvLHrElRqASG+zlg3tQ/k1nX/DiEi49WqZ40plUps2bIFISEh9YYgAAgJCcGPP/4IpVIJtVqNH374AaWlpXj44YfrvaesrAwFBQU1XmS6zKQSLA/3B/C/qdPV7v1425k0PLbuOFJzi5utNtKNWi1g9W+XseSXOKjUAib06YBNM4MZgohMnFEFoUWLFsHW1hbOzs5ITU3Frl27Grx+27ZtqKiogLOzM2QyGZ577jn88ssv8PPzq/ee1atXQy6Xa14eHh76/jbIyIQFKLBuat1Tqb+YGoTvZgfDydYSlzIKMPrTI/gjPlukSqk+pRUqvLD1HL78+zoAYMHwzvjPhJ6wNDeqX4FEZACido0tXrwY7777boPXXL58GV27dgUA3L59G0qlEikpKVixYgXkcjn27t0LiaTu6awvvvgioqOj8c4776Bt27bYuXMnPvzwQxw5cgQ9evSo856ysjKUlf1vZdmCggJ4eHiwa4ygUgv1TqXOzC/Bv7acw7nUPADAcw/5YOEjXWBuxjdaseXeLcPTm88gJjUPlmZSrJ3QE2MD24tdFhEZmLZdY6IGoVu3biE3N7fBa3x8fGBpaVnreFpaGjw8PHD8+HEMGDCg1vmkpCT4+fkhLi4O3bt31xwfNmwY/Pz88MUXX2hVI8cIkbbKK9VY83sCNhxLBgAEezvhs0mBcHXgQGqxJN26i5kbTyNVWQy5tQW+nNYH/X04HZ7IFGj7/i3qOkIuLi5wcXFp1L1qddWO4Pe23tyruLhqrIZUWvNf5GZmZpp7ifTJ0lyKN8L90cfLEa/+fB7RyUqM+uQoPpscyDdfA6urte7MDSWe/e4s8ksq4Olkg40z+8HXxU7sUomohTGKWWOnTp3C6dOnMXDgQDg6OiIpKQnLli1DdnY2Ll26BJlMhvT0dISGhmLz5s0IDg5GRUUF/P39oVAo8N5778HZ2Rk7d+7EwoULsXfvXowaNUqrr80WIWqMpFt38XzkOSRmF0IqARaO6IrnBvtAylWJ9S4qLhMr9sTXWPiyjbUF7pZVolItINCzDb55qi+c7WQiVklEza1VzRqzsbHBjh07EBoaii5dumD27Nno2bMnDh8+DJms6pdbRUUFEhMTNS1BFhYW+O233+Di4oLw8HD07NkTmzdvxrfffqt1CCJqLF8XO/zyrxCMD2wPtQC8G5VQ1TpRXCF2aa1KVFwm5kaeq7X6d15JRVUI8miDrc/0ZwgionoZRYuQmNgiRE0hCAK2Rt/Em7svoVylhoeTNdZN6YOA9vIGB1/T/VWv8dTQFigKuRWOLhrK50pkgoxijBBRayeRSDD5AU/0aC/H3C1ncVNZgvHrjuPxPh1wMCEHWfe8iSvkVlge7s+dzbV0v33gAO4DR0T3ZxRdY0TGrkcHOX59cRBCu7qivFKNLadSa4QgoGoD0LmR5xAVV/du91QT94EjIn1gECJqJnIbC3wxtQ/sZXU3xFb3Ua/YEw+Vmj3WDUm6dRdbT6VqdS33gSOihrBrjKgZnUm5g8KyynrPC2B3TkMy80vw8YGr+Ols2n3DogRVq38Hezs1T3FEZJQYhIiaEbtzGudOUTnWHU7CpuM3UF5ZtQ7YsG6u6O/rjFV7LwNAjU1xq4dGLw/350BpImoQgxBRM9K2m2ZrdCo6udrD3920ZyoWlVVi47FkfHn4uqYlLbijExaN7II+XlUtPR3aWNdaR8iNA8+JSEucPn8fnD5P+lQ95TsrvxTa/MUL7eqK54f4oY+Xo8Fra0nKK9XYGp2KTw9dw+27VavHd1M44NWwLni4s0ut/QW5FAER/ZNR7DVmDBiESN+qFwEE6u7OeX1UN1xIz8evFzJQPQymv48T/jXEDwP92ta7ybCxaCi0qNQCdp9Pxwd/XMFNZQkAwNPJBgse6Yzwnu5cmZuItMYgpCcMQmQIdW0L8c91hJJvF+HLw0nYfi4NFaqqv6Y9O8jx/MN+eMS/nVGGgvq+7zce9YeluRT/2ZeIhKxCAICLvQzzQjthYl8PWJpzgisR6YZBSE8YhMhQtO3Oycwvwdd/J+P76BSUVlQNFO7kaoe5D/siopc7zM2kOn0+sVS3hN3vF469lTnmPuyLGSEdYWPJYYxE1DgMQnrCIEQtRe7dMmw8dgPfnriBwtKqgcMdHK0x5yFfOFiZY/XvCQ22MIlJm+0wAODZwT54/mFftLGxbKbKiKi1YhDSEwYhamkKSisQeTIF648kI7eovN7rqtuC1k0NanQY0lcr04mkXEz6+uR9r9v6TH+un0REesG9xohaKQcrCzz/sB9mhnjjh9OpWLk3HnWtLSigKgyt2BOP4f5uOgcYbcYx1Se/uAKXMvIRl5GPuPQCnLyeq9XX5PpJRNTcGISIjJS1pRm6ujnUGYKqVa9UPfLjv9GjfRt4t7WBl7MtvNvawsvZBvZWFnXeV994nur90O5tZbpVWIa4jHxcSq8KPXEZ+Ui7U9Ko74nbYRBRc2MQIjJi2ragXMm+iyvZd2sdb2tnCS9nW3R0tkVHZxt4tbWFp6MNlu++VOeg5upjr/x0AdtO38SlzAJkF5TV+TU9nWwQ0N4B3d3l8Fc4YNH2C7hVWFbn5+V2GEQkFgYhIiOmbQvKvKF+sDSXIvl2MVJyi3Ajtwi375ZrXmdT7uj0de+WVeJQ4i0AgEQC+LrYobu7AwLc5eje3gHdFXLIbWq2Nr01pjvmRp6DBNwOg4haDgYhIiMW7O0Ehdyq3pWqq1taXhrWuVbIKCytQEpuMW7kFiEltxjJt4uQkluEhMzCBjeGrTYhqAMmPeCBrm4OsJXd/1dJWIAC66YGcTsMImpRGISIjJiZVILl4f6Nammxt7JAQHs5AtrLaxzXdobXY306aPb70lZYgALD/d1a9HpHRGRauFwrkZGrbmlxk9fsJnOTWzVq6nx1K1N90USCqtljjR3PYyaVYICvM8b0bo8Bvs4MQUQkKrYIEbUC+mxpaUorExGRseGCivfBBRXJVDVlHSEiIrFxQUUiahKO5yEiU8AgRET1qh7PQ0TUWnGwNBEREZksBiEiIiIyWQxCREREZLIYhIiIiMhkMQgRERGRyWIQIiIiIpPFIEREREQmi0GIiIiITBaDEBEREZksrix9H9VbsRUUFIhcCREREWmr+n37fluqMgjdR2FhIQDAw8ND5EqIiIhIV4WFhZDL5fWe5+7z96FWq5GRkQF7e3tIJNxs8l4FBQXw8PDAzZs3G9zZl/SLz10cfO7i4HMXR2t47oIgoLCwEO7u7pBK6x8JxBah+5BKpejQoYPYZbRoDg4ORvsXxZjxuYuDz10cfO7iMPbn3lBLUDUOliYiIiKTxSBEREREJotBiBpNJpNh+fLlkMlkYpdiUvjcxcHnLg4+d3GY0nPnYGkiIiIyWWwRIiIiIpPFIEREREQmi0GIiIiITBaDEBEREZksBiG6r/T0dEydOhXOzs6wtrZGjx49cObMGc15QRDwxhtvQKFQwNraGsOGDcPVq1dFrNj4dezYERKJpNbrX//6FwCgtLQU//rXv+Ds7Aw7Ozs89thjyM7OFrlq46dSqbBs2TJ4e3vD2toavr6+WLlyZY29ivjzbhiFhYV4+eWX4eXlBWtra4SEhOD06dOa83zuTff3338jPDwc7u7ukEgk2LlzZ43z2jxjpVKJKVOmwMHBAW3atMHs2bNx9+7dZvwuDEAgaoBSqRS8vLyEGTNmCKdOnRKuX78u7Nu3T7h27ZrmmjVr1ghyuVzYuXOncP78eSEiIkLw9vYWSkpKRKzcuOXk5AiZmZma1x9//CEAEP78809BEARhzpw5goeHh3Dw4EHhzJkzQv/+/YWQkBBxi24FVq1aJTg7Owt79+4VkpOThZ9++kmws7MTPv74Y801/Hk3jCeeeELw9/cXDh8+LFy9elVYvny54ODgIKSlpQmCwOeuD7/99puwZMkSYceOHQIA4ZdffqlxXptnHBYWJvTq1Us4efKkcOTIEcHPz0+YNGlSM38n+sUgRA1atGiRMHDgwHrPq9Vqwc3NTfjPf/6jOZaXlyfIZDJh69atzVGiSXjppZcEX19fQa1WC3l5eYKFhYXw008/ac5fvnxZACCcOHFCxCqN3+jRo4VZs2bVODZ+/HhhypQpgiDw591QiouLBTMzM2Hv3r01jgcFBQlLlizhczeAfwYhbZ5xfHy8AEA4ffq05prff/9dkEgkQnp6erPVrm/sGqMG7d69G3379sXjjz8OV1dXBAYG4uuvv9acT05ORlZWFoYNG6Y5JpfL8cADD+DEiRNilNzqlJeXIzIyErNmzYJEIsHZs2dRUVFR45l37doVnp6efOZNFBISgoMHD+LKlSsAgPPnz+Po0aMYOXIkAP68G0plZSVUKhWsrKxqHLe2tsbRo0f53JuBNs/4xIkTaNOmDfr27au5ZtiwYZBKpTh16lSz16wvDELUoOvXr2PdunXo1KkT9u3bh7lz52LevHn49ttvAQBZWVkAgHbt2tW4r127dppz1DQ7d+5EXl4eZsyYAaDqmVtaWqJNmzY1ruMzb7rFixfjySefRNeuXWFhYYHAwEC8/PLLmDJlCgD+vBuKvb09BgwYgJUrVyIjIwMqlQqRkZE4ceIEMjMz+dybgTbPOCsrC66urjXOm5ubw8nJyaj/P3D3eWqQWq1G37598c477wAAAgMDERcXhy+++ALTp08XuTrTsH79eowcORLu7u5il9Lqbdu2DVu2bMH333+P7t27IzY2Fi+//DLc3d35825g3333HWbNmoX27dvDzMwMQUFBmDRpEs6ePSt2adTKsUWIGqRQKODv71/jWLdu3ZCamgoAcHNzA4BaM5ays7M156jxUlJScODAATz99NOaY25ubigvL0deXl6Na/nMm27hwoWaVqEePXpg2rRpmD9/PlavXg2AP++G5Ovri8OHD+Pu3bu4efMmoqOjUVFRAR8fHz73ZqDNM3Zzc0NOTk6N85WVlVAqlUb9/4FBiBr04IMPIjExscaxK1euwMvLCwDg7e0NNzc3HDx4UHO+oKAAp06dwoABA5q11tZo48aNcHV1xejRozXH+vTpAwsLixrPPDExEampqXzmTVRcXAyptOavRTMzM6jVagD8eW8Otra2UCgUuHPnDvbt24cxY8bwuTcDbZ7xgAEDkJeXV6OV7tChQ1Cr1XjggQeavWa9EXu0NrVs0dHRgrm5ubBq1Srh6tWrwpYtWwQbGxshMjJSc82aNWuENm3aCLt27RIuXLggjBkzhtNa9UClUgmenp7CokWLap2bM2eO4OnpKRw6dEg4c+aMMGDAAGHAgAEiVNm6TJ8+XWjfvr1m+vyOHTuEtm3bCq+++qrmGv68G0ZUVJTw+++/C9evXxf2798v9OrVS3jggQeE8vJyQRD43PWhsLBQiImJEWJiYgQAwgcffCDExMQIKSkpgiBo94zDwsKEwMBA4dSpU8LRo0eFTp06cfo8tX579uwRAgICBJlMJnTt2lX46quvapxXq9XCsmXLhHbt2gkymUwIDQ0VEhMTRaq29di3b58AoM5nWVJSIjz//POCo6OjYGNjI4wbN07IzMwUocrWpaCgQHjppZcET09PwcrKSvDx8RGWLFkilJWVaa7hz7th/Pjjj4KPj49gaWkpuLm5Cf/617+EvLw8zXk+96b7888/BQC1XtOnTxcEQbtnnJubK0yaNEmws7MTHBwchJkzZwqFhYUifDf6IxGEe5ZMJSIiIjIhHCNEREREJotBiIiIiEwWgxARERGZLAYhIiIiMlkMQkRERGSyGISIiIjIZDEIERE1QWJiIvr16wdvb2/s2rVL7HKISEdcR4iIqAkmTpyI4OBg9OzZE7Nnz9bsw0dExoEtQkTUInTs2BEfffSR2GXUKzExEW5ubigsLKxxXC6Xw8vLC35+fnB1dW3wc3z66aeQSCQICQlBcXFxo+r44osvEB4e3qh7iag2BiEiapLw8HCEhYXVee7IkSOQSCS4cOFCM1elf6+99hpefPFF2Nvb1zj+1ltvYeLEifDz88Nrr71W7/1btmzBK6+8gk8++QRKpRKPPfYYKioqalyTm5uLsLAwuLu7QyaTwcPDAy+88AIKCgo018yaNQvnzp3DkSNH9PsNEpkoBiEiapLZs2fjjz/+QFpaWq1zGzduRN++fdGzZ08RKtOf1NRU7N27FzNmzKh17tSpU+jQoQOefPJJHD9+vM77f/vtN8yZMwc//fQTXnzxRfz999/IyMjAjBkzcO/oBKlUijFjxmD37t24cuUKNm3ahAMHDmDOnDmaaywtLTF58mR88sknev8+iUwRgxARNcmjjz4KFxcXbNq0qcbxu3fv4qeffsLs2bMBANu3b0f37t0hk8nQsWNHvP/++/V+zhs3bkAikSA2NlZzLC8vDxKJBH/99RcA4K+//oJEIsG+ffsQGBgIa2trDB06FDk5Ofj999/RrVs3ODg4YPLkyTW6odRqNVavXg1vb29YW1ujV69e+Pnnnxv8Hrdt24ZevXqhffv2tc5t3LgRkydPxrRp0xAZGYnKysoa548dO4bp06djx44diIiIAAC4urrir7/+wrVr1/DSSy9prnV0dMTcuXPRt29feHl5ITQ0FM8//3yt1p/w8HDs3r0bJSUlDdZNRFoQdctXImoVFi5cKPj6+gpqtVpzbMOGDYK1tbWQl5cnnDlzRpBKpcJbb70lJCYmChs3bhSsra2FjRs3aq738vISPvzwQ0EQBCE5OVkAIMTExGjO37lzRwAg/Pnnn4Ig/G8n7f79+wtHjx4Vzp07J/j5+QkPPfSQ8Mgjjwjnzp0T/v77b8HZ2VlYs2aN5vO8/fbbQteuXYWoqCghKSlJ2LhxoyCTyYS//vqr3u8vIiJCmDNnTq3j2dnZgoWFhRAXFydUVlYK7dq1E3bu3Nm4h1iH9PR04aGHHhKmTJlS43hRUZEglUo1z4KIGo8tQkTUZLNmzUJSUhIOHz6sObZx40Y89thjkMvl+OCDDxAaGoply5ahc+fOmDFjBl544QX85z//afLXfvvtt/Hggw8iMDAQs2fPxuHDh7Fu3ToEBgZi0KBBmDBhAv78808AQFlZGd555x1s2LABI0aMgI+PD2bMmIGpU6fiyy+/rPdrpKSkwN3dvdbxyMhIdO/eHd27d4eZmRmefPLJWi1jjTFp0iTY2Nigffv2cHBwwDfffFPjvI2NDeRyOVJSUpr8tYhMHYMQETVZ165dERISgg0bNgAArl27hiNHjmi6xS5fvowHH3ywxj0PPvggrl69CpVK1aSvfe/4o3bt2sHGxgY+Pj41juXk5GjqKi4uxvDhw2FnZ6d5bd68GUlJSfV+jZKSElhZWdU6vnHjRkydOlXz8dSpU/Hrr7/i1q1bTfqePvzwQ5w7dw67du1CUlIS/v3vf9e6xtrautEzz4jof8zFLoCIWofZs2fjxRdfxOeff46NGzfC19cXDz30UKM+l1Ra9W804Z6BxP+cYVXNwsJC82eJRFLj4+pjarUaQNW4JQD49ddfa433kclk9dbTtm1b3Llzp8axM2fOIC4uDq+++ioWLVqkOa5SqRAZGYn58+fX+/nux83NDW5ubujatSucnJwwaNAgLFu2DAqFQnONUqmEi4tLo78GEVVhixAR6cUTTzwBqVSK77//Hps3b8asWbMgkUgAAN26dcOxY8dqXH/s2DF07twZZmZmtT5X9Rt8Zmam5ti9A6cby9/fHzKZDKmpqfDz86vx8vDwqPe+wMBAxMfH1zi2ceNGDB48GOfPn0dsbKzm9eqrr+qle6xadYgrKyvTHEtKSkJpaSkCAwP19nWITBVbhIhIL+zs7DBx4kS89tprKCgoqDHVfMGCBejXrx9WrlyJiRMn4sSJE/jss8/w3//+t87PZW1tjf79+2PNmjXw9vZGTk4Oli5d2uQa7e3t8corr2D+/PlQq9UYOHAg8vPzcezYMTg4OGD69Ol13jdixAg8/fTTUKlUMDMzQ1lZGbZu3Yp33nkHAQEBNa59+umnsXbtWpw7dw5BQUE61ffbb78hOzsb/fr1g52dHS5duoSFCxfiwQcfRMeOHTXXHTlyBD4+PvD19dX5GRBRTWwRIiK9mT17Nu7cuYMRI0bUGFwcFBSEbdu24YcffkBAQADeeOMNvPXWW3Wuy1Ntw4YNqKysRJ8+ffDyyy/j7bff1kuNK1euxLJly7B69Wp069YNYWFh+PXXX+Ht7V3vPSNHjoS5uTkOHDgAANi5cyfy8/Mxbty4Wtd26tQJPXr0wMaNG3WuzdraGl9//TUGDhyIbt26Yf78+YiIiMDevXtrXLd161Y888wzOn9+IqqNe40REWnh888/x+7du7Fv3z5R67h06RKGDh2KK1euQC6Xi1oLUWvArjEiIi0899xzyMvLQ2FhYa1tNppTZmYmNm/ezBBEpCdsESIiIiKTxTFCREREZLIYhIiIiMhkMQgRERGRyWIQIiIiIpPFIEREREQmi0GIiIiITBaDEBEREZksBiEiIiIyWQxCREREZLIYhIiIiMhk/R86NXVuSU9qEQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "energies = []\n",
    "volumes = []\n",
    "\n",
    "structure_copy = structure.copy()\n",
    "original_cell = structure.cell.copy()\n",
    "\n",
    "structure_copy.calc = calc\n",
    "\n",
    "for scale in np.arange(0.9, 1.11, 0.01):\n",
    "    structure_copy.set_cell(scale * original_cell, scale_atoms=True)\n",
    "    volumes.append(structure_copy.get_volume())\n",
    "    energies.append(structure_copy.get_potential_energy() / len(structure_copy))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(volumes, energies, '-o')\n",
    "ax.set_xlabel('Volume (Å^3)')\n",
    "ax.set_ylabel('Energy (eV/atom)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18e8a3c7",
   "metadata": {},
   "source": [
    "### Calculate dipole moments\n",
    "\n",
    "The calculator also supports the calculation of dipole moments. Note that you must have a NEP model trained to predict dipoles; otherwise the output will be nonsensical. Refer to the [GPUMD documentation](https://gpumd.org/) for information on how to train a dipole model. \n",
    "\n",
    "The model and structure are included in the `calorine/tutorials/` folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7f742e44-de81-4027-a5b7-f6c0110cbe33",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DFT dipole moment:            [4.85503 0.07166 0.16004] a.u.\n",
      "NEP predicted dipole moment:  [4.78034782 0.03274603 0.04843106] a.u.\n"
     ]
    }
   ],
   "source": [
    "dipole_model = 'nep-PTAF-dipole.txt'\n",
    "structure = read('PTAF-monomer.xyz')\n",
    "calc = CPUNEP(dipole_model)\n",
    "structure.calc = calc\n",
    "\n",
    "dipole = structure.get_dipole_moment()\n",
    "dft_dipole = structure.info['dipole']\n",
    "\n",
    "print('DFT dipole moment:'.ljust(30) +  f'{dft_dipole} a.u.')\n",
    "print(f'NEP predicted dipole moment:'.ljust(30) + f'{dipole} a.u.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5ed3304-a3b0-4ed8-9e1d-166355ab2a74",
   "metadata": {
    "tags": []
   },
   "source": [
    "## GPU-based calculator\n",
    "\n",
    "The basic usage of the `GPUNEP` calculator class is completely analoguous to the `CPUNEP` calculator.\n",
    "We create an atomic structure and attach the calculator object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6e0bc508-b926-4312-a3ef-43a2f6733be4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from calorine.calculators import GPUNEP\n",
    "\n",
    "structure = bulk('PbTe', crystalstructure='rocksalt', a=6.7)\n",
    "calc = GPUNEP('nep-PbTe.txt')\n",
    "structure.calc = calc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0daccfc5-91ce-4295-b3b5-eb7ab2a1be7c",
   "metadata": {},
   "source": [
    "Afterwards we can readily obtain energies, forces, and stresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "86a7f3fb-3d7f-4b6c-b47c-d4203ad8a739",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy (eV): -7.6631630659\n",
      "Forces (eV/Å):\n",
      " [[ 2.91459e-07  2.77747e-07  1.82196e-07]\n",
      " [-2.91459e-07 -2.77747e-07 -1.82196e-07]]\n",
      "Stresses (eV/Å^3):\n",
      " [ 1.09085639e-02  1.09085643e-02  1.09085477e-02 -3.55200559e-10\n",
      " -3.55200534e-10  9.51932078e-10]\n"
     ]
    }
   ],
   "source": [
    "print('Energy (eV):', structure.get_potential_energy())\n",
    "print('Forces (eV/Å):\\n', structure.get_forces())\n",
    "print('Stresses (eV/Å^3):\\n', structure.get_stress())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec80392-5346-4da3-9cf1-c556f31a2e37",
   "metadata": {},
   "source": [
    "### Temporary and specified directories\n",
    "\n",
    "Under the hood, the `GPUNEP` calculator creates a directory and writes the input files necessary to run GPUMD.\n",
    "By default, this is done in temporary directories that are automatically removed once the calculations has finished.\n",
    "It is also possible to run in a user-specified directory that will be kept after the calculations finish.\n",
    "This is especially useful when running molecular dynamics simulations (see below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be3a38a4-7d18-43ec-b512-37ce7a973223",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.6631630659"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc.set_directory('my_directory')\n",
    "structure.get_potential_energy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eaa6432-df7c-4de5-bb66-f299bf002c6a",
   "metadata": {},
   "source": [
    "After this is run, there should be a new directory, `my_directory`, in which input and output files from GPUMD are available.\n",
    "This can be useful for debugging."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0003105-86a3-4253-94c5-26f10c27dbc7",
   "metadata": {},
   "source": [
    "### Running custom molecular dynamics simulations\n",
    "To take advantage of the Python workflow as well as raw speed of the GPU accelerated NEP implementation, the `GPUNEP` calculator contains a convenience function for running customized molecular dynamics simulations.\n",
    "This should typically be done in a specified directory.\n",
    "The [parameters of the run.in file](https://gpumd.zheyongfan.org/index.php/Main_Page#Inputs_for_the_src.2Fgpumd_executable) are specified as a list of tuples with two elements, the first being the keyword name and the second any arguments to that keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f8e8954f-cebb-4d22-a0b6-670603d03ce9",
   "metadata": {},
   "outputs": [],
   "source": [
    "calc.set_directory('my_md_simulation')\n",
    "supercell = structure.repeat(5)\n",
    "structure.calc = calc\n",
    "parameters = [('velocity', 300),\n",
    "              ('dump_thermo', 100),\n",
    "              ('dump_position', 100),\n",
    "              ('ensemble', ('nvt_lan', 300, 300, 100)),\n",
    "              ('run', 1000)]\n",
    "calc.run_custom_md(parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6b6f5e1-c44c-45d1-bb5d-f522393a8cd2",
   "metadata": {},
   "source": [
    "Once this is run, the results are available in the folder `my_md_simulations`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "6be099dd3823815779ed5b2d986ab000320f08fb94733d934e13608f08061690"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
